import os
import glob
import pandas as pd
from matplotlib.gridspec import GridSpec
# Import cerebellum packages
import matplotlib.pyplot as plt
import SUITPy.flatmap as flatmap
bids_dir = '/data/projects/social_doors/'
os.chdir(bids_dir)
data_dir = os.path.join(bids_dir, 'derivatives','social_doors-nilearn')
# Define subject list
#subjs_scan_info = pd.read_csv(bids_dir+'/derivatives/mriqc/mriqc_summary_poor.csv')
#subjs_list = list(subjs_scan_info['subject'].unique())
#subjs_list.sort()
subjs_info = pd.read_csv(bids_dir+'derivatives/participants_good.tsv', sep='\t', index_col=0)
subjs_info = subjs_info.rename(columns={'participant_id': 'subject_label'})
# Remove participants with bad data
subjs_info = subjs_info[subjs_info['subject_label'].str.contains('sub-3880')==False]
subjs_info = subjs_info[subjs_info['subject_label'].str.contains('sub-4069')==False]
subjs_info_kids = subjs_info[subjs_info['group']=='kid']
subjs_info_colg = subjs_info[subjs_info['group']=='college']
subjs_list = subjs_info['subject_label'].to_list()
subjs_list_kids = subjs_info_kids['subject_label'].to_list()
subjs_list_colg = subjs_info_colg['subject_label'].to_list()
print('Found '+str(len(subjs_list_kids))+' adolescent subjects')
print('Found '+str(len(subjs_list_colg))+' college subjects')
Found 32 adolescent subjects Found 29 college subjects
# Filter for releveant data
subjs_info_fltr = subjs_info[['subject_label', 'age']]
subjs_info_kids_fltr = subjs_info_kids[['subject_label','age']]
subjs_info_colg_fltr = subjs_info_colg[['subject_label','age']]
from nilearn.glm.second_level import make_second_level_design_matrix
design_matrix = make_second_level_design_matrix(subjs_list_kids, subjs_info_kids_fltr)
/opt/anaconda3/lib/python3.7/site-packages/nilearn/glm/__init__.py:56: FutureWarning: The nilearn.glm module is experimental. It may change in any future release of Nilearn.
'It may change in any future release of Nilearn.', FutureWarning)
/opt/anaconda3/lib/python3.7/site-packages/nilearn/glm/first_level/design_matrix.py:475: UserWarning: Attention: Design matrix is singular. Aberrant estimates are expected.
warn('Attention: Design matrix is singular. Aberrant estimates '
from nilearn.plotting import plot_design_matrix
plot_design_matrix(design_matrix)
<AxesSubplot:label='conditions', ylabel='scan number'>
from nilearn.glm.second_level import SecondLevelModel
mni_mask = bids_dir+"derivatives/fmriprep/sub-010/anat/sub-010_space-MNI152NLin2009cAsym_label-GM_probseg_bin.nii.gz"
contrast = 'valence_x_outcome'
task = 'social'
group = 'kids'
temp_file_list = []
for subj in subjs_list_kids:
temp_file = glob.glob(os.path.join(data_dir,subj,'zmap_'+task+'_'+contrast+'.nii'))
temp_file_list.append(temp_file[0])
model = SecondLevelModel(mask_img=mni_mask, smoothing_fwhm=6.0)
model.fit(temp_file_list, design_matrix=design_matrix)
SecondLevelModel(mask_img='/data/projects/social_doors/derivatives/fmriprep/sub-010/anat/sub-010_space-MNI152NLin2009cAsym_label-GM_probseg_bin.nii.gz',
smoothing_fwhm=6.0)
z_map = model.compute_contrast('intercept', output_type='z_score')
z_map.to_filename(os.path.join(data_dir,'group',
'zmap_'+group+'_'+task+'_'+contrast+'.nii.gz'))
from nilearn.plotting import plot_stat_map
plot_stat_map(z_map, threshold=2, cut_coords=range(-65,66,10), display_mode='x')
/opt/anaconda3/lib/python3.7/site-packages/nilearn/plotting/img_plotting.py:300: FutureWarning: Default resolution of the MNI template will change from 2mm to 1mm in version 0.10.0 anat_img = load_mni152_template()
<nilearn.plotting.displays._slicers.XSlicer at 0x7f1b4c0b8390>
from nilearn.glm import threshold_stats_img
_, threshold = threshold_stats_img(z_map, alpha=0.001, height_control='fpr')
plot_stat_map(z_map, threshold=threshold, cut_coords=range(-65,66,10), display_mode='x')
<nilearn.plotting.displays._slicers.XSlicer at 0x7f1b40b60390>
design_matrix = make_second_level_design_matrix(subjs_list_kids, subjs_info_kids_fltr)
#contrasts = ['valence_x_outcome',
# 'positive','positive_win','positive_loss',
# 'negative','negative_win','negative_loss']
contrasts=['all_winVlos', 'positive_winVlos']
contrasts_2nd = ['age','intercept']
tasks = ['mdoors','social']
group = 'kids'
alpha = 0.001
threshold_data = {}
for task in tasks:
for contrast in contrasts:
temp_file_list = []
for subj in subjs_list_kids:
temp_file = glob.glob(os.path.join(data_dir,subj,'zmap_'+task+'_'+contrast+'.nii'))
temp_file_list.append(temp_file[0])
temp_file_list.sort()
print('Calculating '+task+' group '+contrast+' contrast')
model = SecondLevelModel(mask_img=mni_mask, smoothing_fwhm=8.0)
model.fit(temp_file_list, design_matrix=design_matrix)
for contrast_2nd in contrasts_2nd:
z_map = model.compute_contrast(contrast_2nd, output_type='z_score')
z_map.to_filename(os.path.join(data_dir,'group',
'zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_unc.nii.gz'))
# Multiple Comparisons Correction
z_map_thresh, threshold = threshold_stats_img(z_map, alpha=alpha, height_control='fpr')
z_map_thresh.to_filename(os.path.join(data_dir,'group',
'zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_fpr.nii.gz'))
threshold_data['zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_fpr_'+str(alpha)] = threshold
Calculating mdoors group all_winVlos contrast Calculating mdoors group positive_winVlos contrast Calculating social group all_winVlos contrast Calculating social group positive_winVlos contrast
def plot_stat_map_cb(filename, threshold):
# Find stat map
stat_filename = os.path.join(data_dir,'group',filename+'.nii.gz')
# Import stat map as a cerebellum flatmap
funcdata = flatmap.vol_to_surf(stat_filename,
space='SUIT')
# Set figure specs
fig = plt.figure(figsize=(15, 4))
gs = GridSpec(2, 3)
ax_img1 = plt.subplot(gs[0, :2])
ax_img2 = plt.subplot(gs[1, :2])
# Plot stat maps
plot_stat_map(stat_filename, threshold=threshold, cut_coords=range(-65,0,10), display_mode='x', axes=ax_img1,
annotate=False, title='Left hemisphere')
plot_stat_map(stat_filename, threshold=threshold, cut_coords=range(66,5,-10), display_mode='x', axes=ax_img2,
annotate=False, title='Right hemisphere')
ax_joint = plt.subplot(gs[:, 2:])
ax_joint.set(title='Cerebellum flatmap')
flatmap.plot(data=funcdata, cmap='hot',
threshold=[threshold, threshold+5],
colorbar=True,
render='matplotlib', new_figure=False)
filename = 'zmap_kids_mdoors_positive_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]
plot_stat_map_cb(filename+'unc', threshold)
filename = 'zmap_kids_mdoors_all_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]
plot_stat_map_cb(filename+'unc', threshold)
filename = 'zmap_kids_social_positive_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]
plot_stat_map_cb(filename+'unc', threshold)
filename = 'zmap_kids_social_all_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]
plot_stat_map_cb(filename+'unc', threshold)
design_matrix = make_second_level_design_matrix(subjs_list_colg, subjs_info_colg_fltr)
contrasts=['all_winVlos', 'positive_winVlos']
contrasts_2nd = ['age','intercept']
tasks = ['mdoors','social']
group = 'colg'
alpha = 0.001
threshold_data = {}
for task in tasks:
for contrast in contrasts:
temp_file_list = []
for subj in subjs_list_colg:
temp_file = glob.glob(os.path.join(data_dir,subj,'zmap_'+task+'_'+contrast+'.nii'))
temp_file_list.append(temp_file[0])
temp_file_list.sort()
print('Calculating '+task+' group '+contrast+' contrast')
model = SecondLevelModel(mask_img=mni_mask, smoothing_fwhm=8.0)
model.fit(temp_file_list, design_matrix=design_matrix)
for contrast_2nd in contrasts_2nd:
z_map = model.compute_contrast(contrast_2nd, output_type='z_score')
z_map.to_filename(os.path.join(data_dir,'group',
'zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_unc.nii.gz'))
# Multiple Comparisons Correction
z_map_thresh, threshold = threshold_stats_img(z_map, alpha=alpha, height_control='fpr')
z_map_thresh.to_filename(os.path.join(data_dir,'group',
'zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_fpr.nii.gz'))
threshold_data['zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_fpr_'+str(alpha)] = threshold
/opt/anaconda3/lib/python3.7/site-packages/nilearn/glm/first_level/design_matrix.py:475: UserWarning: Attention: Design matrix is singular. Aberrant estimates are expected.
warn('Attention: Design matrix is singular. Aberrant estimates '
Calculating mdoors group all_winVlos contrast Calculating mdoors group positive_winVlos contrast Calculating social group all_winVlos contrast Calculating social group positive_winVlos contrast
filename = 'zmap_colg_mdoors_positive_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]
plot_stat_map_cb(filename+'unc', threshold)
filename = 'zmap_colg_mdoors_all_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]
plot_stat_map_cb(filename+'unc', threshold)
filename = 'zmap_colg_social_positive_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]
plot_stat_map_cb(filename+'unc', threshold)
filename = 'zmap_colg_social_all_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]
plot_stat_map_cb(filename+'unc', threshold)
subjs_info_num = subjs_info.copy()
subjs_info_num = subjs_info_num.replace({'sex': {'F': 0, 'M': 1},
'group': {'college': -1, 'kid': 1}})
subjs_info_num = subjs_info_num.drop(columns=['age', 'sex'])
design_matrix = make_second_level_design_matrix(subjs_list, subjs_info_num)
plot_design_matrix(design_matrix)
<AxesSubplot:label='conditions', ylabel='scan number'>
#design_matrix = make_second_level_design_matrix(subjs_list, subjs_info_fltr)
#contrasts = ['valence_x_outcome',
# 'positive','positive_win','positive_loss',
# 'negative','negative_win','negative_loss']
contrasts=['all_winVlos', 'positive_winVlos']
contrasts_2nd = ['intercept','group']
tasks = ['mdoors','social']
group = 'all'
alpha = 0.001
threshold_all = {}
#temp_file_list = []
#for subj in subjs_list_kids:
# temp_file = glob.glob(os.path.join(data_dir,subj,'zmap_'+task+'_'+contrast+'.nii.gz'))
# temp_file_list.append(temp_file[0])
for task in tasks:
for contrast in contrasts:
temp_file_list = []
for subj in subjs_list:
temp_file = glob.glob(os.path.join(data_dir,subj,'zmap_'+task+'_'+contrast+'.nii'))
temp_file_list.append(temp_file[0])
temp_file_list.sort()
print('Calculating '+task+' group '+contrast+' contrast')
model = SecondLevelModel(mask_img=mni_mask, smoothing_fwhm=8.0)
model.fit(temp_file_list, design_matrix=design_matrix)
for contrast_2nd in contrasts_2nd:
z_map = model.compute_contrast(contrast_2nd, output_type='z_score')
z_map.to_filename(os.path.join(data_dir,'group',
'zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_unc.nii.gz'))
# Multiple Comparisons Correction
z_map_thresh, threshold = threshold_stats_img(z_map, alpha=alpha, height_control='fpr')
z_map_thresh.to_filename(os.path.join(data_dir,'group',
'zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_fpr.nii.gz'))
threshold_data['zmap_'+group+'_'+task+'_'+contrast+'_'+contrast_2nd+'_fpr_'+str(alpha)] = threshold
Calculating mdoors group all_winVlos contrast Calculating mdoors group positive_winVlos contrast Calculating social group all_winVlos contrast Calculating social group positive_winVlos contrast
Where adolescents > young adults
filename = 'zmap_all_social_positive_winVlos_group_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]
plot_stat_map_cb(filename+'unc', threshold)
filename = 'zmap_all_social_all_winVlos_group_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]
plot_stat_map_cb(filename+'unc', threshold)
filename = 'zmap_all_mdoors_positive_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]
plot_stat_map_cb(filename+'unc', threshold)
filename = 'zmap_all_mdoors_positive_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]
plot_stat_map_cb(filename+'unc', threshold)
filename = 'zmap_all_social_positive_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]
plot_stat_map_cb(filename+'unc', threshold)
filename = 'zmap_all_mdoors_all_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]
plot_stat_map_cb(filename+'unc', threshold)
filename = 'zmap_all_social_all_winVlos_intercept_'
threshold = threshold_data[filename+'fpr_'+str(alpha)]
plot_stat_map_cb(filename+'unc', threshold)
behav_data = pd.read_excel(bids_dir+'derivatives/behavioral/aggregate_full_data_SPSS_020419.xlsx',
sheet_name='Sheet1')
behav_data.head()
| Order | Subject | ch_totanx | ch_socialanx | p_totanx | p_socialanx | bfne | ch_cdi_total | ch_cdi_interpersonal | ch_cdi_anhedonia | ... | FB_VS_6mm_R_Like_Lose_NZMean_1 | FB_VS_6mm_R_Like_Lose_NZCount_1 | FB_VS_6mm_R_Like_Lose_NZMin_1 | FB_VS_6mm_R_Like_Lose_NZMax_1 | FB_VS_6mm_R_Like_Win_Mean_1 | FB_VS_6mm_R_Like_Win_NZMean_1 | FB_VS_6mm_R_Like_Win_NZCount_1 | FB_VS_6mm_R_Like_Win_NZMin_1 | FB_VS_6mm_R_Like_Win_NZMax_1 | FB_VS6mm_L_Correct_DislikeVsLike | |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0 | 1 | s010 | 51 | 14 | 7 | 1 | 33 | 11.0 | 0.0 | 5.0 | ... | 0.029036 | 123.0 | -0.296210 | 0.647286 | 0.256033 | 0.256033 | 123.0 | -0.022722 | 0.752379 | -0.127002 |
| 1 | 1 | s011 | 12 | 3 | 4 | 1 | 31 | 8.0 | 0.0 | 3.0 | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
| 2 | 1 | s013 | 8 | 4 | 1 | 1 | 13 | 5.0 | 0.0 | 2.0 | ... | -0.308424 | 123.0 | -0.554069 | 0.030482 | 0.093248 | 0.093248 | 123.0 | -0.183241 | 1.003483 | -0.167305 |
| 3 | 1 | s015 | 18 | 5 | 1 | 0 | 21 | NaN | NaN | NaN | ... | -0.798996 | 123.0 | -1.841710 | -0.031699 | 0.004660 | 0.004660 | 123.0 | -1.020790 | 0.809193 | 0.485397 |
| 4 | 2 | s024 | 21 | 10 | 0 | 0 | 28 | 1.0 | 0.0 | 0.0 | ... | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
5 rows × 1551 columns
behav_data_fltr = behav_data[['Subject','ch_totanx', 'ch_cdi_total']]
behav_data_fltr['Subject'] = behav_data_fltr['Subject'].str.replace('s','sub-')
behav_data_fltr = behav_data_fltr.rename(columns={'Subject':'subject_label'})
behav_data_fltr.head()
/opt/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:3: SettingWithCopyWarning: A value is trying to be set on a copy of a slice from a DataFrame. Try using .loc[row_indexer,col_indexer] = value instead See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy This is separate from the ipykernel package so we can avoid doing imports until
| subject_label | ch_totanx | ch_cdi_total | |
|---|---|---|---|
| 0 | sub-010 | 51 | 11.0 |
| 1 | sub-011 | 12 | 8.0 |
| 2 | sub-013 | 8 | 5.0 |
| 3 | sub-015 | 18 | NaN |
| 4 | sub-024 | 21 | 1.0 |
# Add behavioral data to design matrix
subjs_info_kids_behav = subjs_info_kids_fltr.merge(behav_data_fltr, on='subject_label')
design_matrix = make_second_level_design_matrix(subjs_list_kids, subjs_info_kids_behav)
plot_design_matrix(design_matrix)
<AxesSubplot:label='conditions', ylabel='scan number'>